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ABSTRACT 

The role of radial velocity (RV) jitter in extrasolar planet search surveys is discussed. 
Based on the maximum likelihood principle, improved statistical algorithms for RV 
fitting and period search are developed. These algorithms incorporate a built-in jitter 
determination, so that resulting estimations of planetary parameters account for this 
jitter automatically. This approach is applied to RV data for several extrasolar plan- 
etary systems. It is shown that many RV planet search surveys suffer from periodic 
systematic errors which increase effective RV jitter and can lead to erroneous conclu- 
sions. For instance, the planet candidate HD74156 d may be a false detection made 
due to annual systematic errors. 

Key words: methods: data analysis - methods: statistical - surveys - techniques: 
radial velocities - stars: planetary systems - stars: individual: HD74156 



1 INTRODUCTION 

When analysing radial velocity (RV) data from planet search 
surveys, we should bear in mind that total uncertainties 
, of these RV measurements are assembled from instrumen- 
' tal uncertainties and a 'jitter'. Partly, this jitter is pro- 
\ duced by various processes on the star leading to instabili- 
. ties of the observed radial velocity. Estimations of planetary 

■ masses and orbital elements depend on full RV uncertain- 
, ties, hence RV jitter should be accounted for in the data 

■ analysis. Usually, empirical models based on a set of stel- 
' lar characteristics are used to asse ss RV jitter (e.g.. IWrighd 
' I2OO5I : ISaar. Butler fc Marcvlll99"3 ). Unfortunately, this way 

of jitter estimation allows accuracies of ~ 50% only or even 
worse. Often, the RV jitter remains almost unconstrained a 
priori (in comparison with instrumental errors) and repre- 
sents an extra unknown parameter. 

It is worth stressing that the jitter also depends on the 
instrument, on the way of observations and obtaining final 
RV measurements. For instance, sufficiently long exposures 
average out stellar oscillations and decrease the apparent RV 
jitter. Extra systematic errors (which have not yet been in- 
vestigated in detail in planet search surveys) should increase 
it. When performing a joint analysis of data from different 
observatories, we must not forget that their effective RV 
jitter may be quite different, implying different statistical 
weights to their RV data. 

Until accurate a priori estimations of RV jitter are con- 
structed, we need to use some statistical algorithm of data 
analysis, which could account properly for the presence of 
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poorly known RV jitter. It is possible to construct a sta- 
tistical jitter estimation based on the scattering of the data 
around the RV model for a given star. The aim of this paper 
is to propose efficient tools implementing this idea. The al- 
gorithm can be organised so that the jitter estimation is au- 
tomatically accounted for in estimations of planetary masses 
and orbital parameters (and vice versa). 

In Section [2] the background connected with RV jitter 
is outlined. In Sections |3] and 3] the traditionally used al- 
gorithms of RV curve fitting and posterior empirical jitter 
determination are briefiy discussed and are shown to be un- 
suitable for our goals. In Section [SI the maximum likelihood 
approach is proposed for joint estimation of RV jitter and pa- 
rameters of the RV curve. It is shown that this approach can 
take into account the presence of unknown RV jitter prop- 
erly. In Section[Sl a modification of the likelihood function is 
introduced. This modification allows to perform a 'preven- 
tive' reduction of the statistical bias in the RV jitter. Several 
other issues connected with biasing of estimations are also 
discussed in this section. In Section [T] the effect of possible 
non-Gaussian distribution of RV errors is considered. It is 
shown that many important properties of the maximum like- 
lihood algorithm constructed in the paper, are not destroyed 
by non-Gaussian nature of RV errors. In Section [S] an effi- 
cient numerical implementation of the analytic algorithm 
is described. This implementation is based on the common 
non-finear Levenberg-Marquardt-Gauss least squares algo- 
rithm. In Section [9] the modified likelihood ratio test is 
proposed for checking consistency of RV models emerging 
in planet searches with RV data. This test incorporates a 
built-in estimation of the RV jitter. Based on this test, a 
generalisation of the Lomb-Scargle periodogram is proposed 
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in Section [To] Results obtained for several RV datasets from 
current planet search surveys are presented in Section [TT] 



2 RADIAL VELOCITY JITTER 

High-precision RV data from planet search surveys suffer 
from the phenomenon called 'RV jitter'. RV measurements 
often shows scattering far beyond the level which is expected 
from their internal uncertainties. Let Vi denote RV mea- 
surements made at the epochs ti. Denoting the internal stan- 
dard errors of Vi as (Tmeas.i, the total variances of RV error 
are usually derived as 

f^i = O-moas.i + O-J, (1) 

where the constant term ctJ, softening differences between 
ai, characterizes the RV jitter. 

It is necessary to clarify the notion 'jitter'. We will name 
the term crj in ((T} as 'jitter' (or 'RV jitter', 'fuU RV jitter') 
regardless its physical nature. In the astrophysical part, the 
RV jitter is inspired by various processes in the star leading 
to an apparent instability of its radial velocity. Also, it de- 
pends on the instrument, on the way of observation and its 
reduction to final radial velocity measurement. For example, 
an exposure as long as 20 — 30 min averages out apparent RV 
variations inspired by stellar oscillations , which have peri- 
ods of several minutes for solar- like stars (|Mavor et al ] |2003l : 
lO'Toole. Tinnev fc Joii3 120081 ) . This decreases the astro- 
physical part of the full RV jitter. However, the astrophysi- 
cal jitter does not represent the only source of RV variations 
beyond the expected noise level. Other sources like extra 
systematic RV errors lie in the instrumentation and in the 
data reduction (but they may depend on stellar properties 
as well). In Section[TT]we will see that effective RV jitter may 
be quite different for different observatories. Note that im- 
perfection of RV models (say, extra Doppler variability due 
to undetected planets in the system) also increase the full 
jitter, but this increase does not depend on an instrument. 

It is worth stressing that we are not intending to find 
here any temporal RV model for the jitter. We model the 
RV jitter in the statistical sense, using the square-additive 
model of RV uncertainties. This simplification should yield 
reliable results for the case, when the jitter has roughly uni- 
form frequency spectrum in the frequency range that we 
are interested in. In planet searches, we are mostly inter- 
ested in periods of RV variations from days to years. This 
m eans that, for examp le, t he stellar oscillations investigated 
bv lMavor et"al] (|2003l ) and lO'Toole. Tinnev fc Joned (|2008l ) 
quite can be processed in this way, because the range of their 
periods (minutes or even hours) lies far beyond the period 
range that we deal with. However, some kinds of extra RV 
variability in the data may require an explicit representa- 
tion in the temporal model of the RV curve. These include, 
for example, quasi-periodic long-period instrumental errors 
(see Section [TTI). RV drifts inspired by spots on the rotati ng 
stellar surface l|Bonfils et al.ll2007l : [Saar fc Donahuelll997l '). 

3 LEAST SQUARES APPROACH 

If we knew the exact statistical weights of the observations, 
Wi oc cr~^, we could write down the full variances of Vi as 



(7i = hi/Wi, (2) 

where the parameter k (the error variance for the unit 
weight) is unspecified. This is the framework which is typ- 
ically assumed in usual statistical algorithms. Clearly, the 
models ((2]) and ([1]) are different. Hereafter, we will refer 
to IT]) as to the 'square-additive' model and to ([2} as to the 
'multiplicative' one. Both these models require a sequence 
of A'' a priori fixed quantities ((Tmcas.i or Wi) and contain an 
unknown 'variance' parameter (cr* or k). If all instrumental 
uncertainties are equal to each other then the models ((TJ 
and ([2]) become equivalent. 

We need to fit our RV observations by a model 
V — ij,{t, 0) which depends on d free parameters forming 
the vector 0. Traditionally, the best-fitting estimations 9* 
are obtained in result of minimizing the function — 
((f — ^)^/(T^) by 00 This is equivalent to minimizing the 
function = Ky^ — (w{v — fj,)^) which does not contain 
any undefined quantities. This is the essence of the least 
squares principle commonly used to obtain the best-fitting 
values of unknown parameters of the RV curve. 

The least squares approach assumes that the weights of 
observations and, hence, the RV jitter are known a priori. 
This a priori ji tter estimation is usually obtained from e m- 
pirical models (|Saar. Butler fc MarcvUToOSl : IWrightll2005l ) or 
even is neglected. Inaccurate values of the jitter inject extra 
bieis in the least squares estimations and decrease their reli- 
ability, especially for the cases when the planetary orbits are 
not constrained well. Still, the accuracy of the a priori jitter 
estimations is not better than ~ 50%. The jitter of several 
m/s (that is, of the order of typical internal RV precision 
reached in planet search surveys) have the largest effect on 
the best-fitting parameters of the RV curve. Unfortunately, 
it is the region where the a priori RV jitter estimations are 
mostly uncertain. 



4 METHOD-OF-MOMENTS ESTIMATOR 

If the true values of the parameters 9 of the RV curve were 
somehow known, we could estimate RV jitter based on the 
observed scattering of the residuals around the RV model 
IJ,{t, 9) as follows: 

= {{v- ^{t, 9)f) /N - (aLas) /N. (3) 

It easy to see that the first term in the right hand side of this 
equation represents the second sample moment of the resid- 
uals. Its mathematical expectation is crj-|-((7^cas)/-^- There- 
fore, the estimator Q could be obtained after equating the 
second sample moment to its expectation. Such an estima- 
tor is called the method-of- moments estimator (MME) . The 
jitter estimation is, probably, the most easy and intu- 
itive one. However, it does not estimate anything but the 
RV jitter. Eventually, we are interested in the estimation of 
9, taking into account some most suitable value of the RV 
jitter. It is possible to organise an iterative process based on 
the MME of the RV jitter and on the least squares estima- 
tor of 9, but this way requires intensive calculations due to 



^ See Appendix |A] for explanations of several mathematical no- 
tations (like the operation (*)) used in the paper. 
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multiple non-linear minimizations and thus is not prac- 
tical. In addition, estimators constructed using the method 
of moments do not necessary provide the best accuracy. It 
is not hard to show that the variance of the MME Q is 
2{a'^)/N'^. As will be shown in Section (5] this is not the 
minimum variance possible for estimating the RV jitter. 



5 MAXIMUM LIKELIHOOD ESTIMATOR 

We note that the least squares principle is often considered 
as a special case of a more general maximum likelihood prin- 
ciple. Assuming that the errors of the RV measurements are 
uncorrelated and Gaussian, we can write down the associ- 
ated log-likelihood function as 

InC^ (Ina) + N\nV2^. (4) 

This function depends on the parameters and k. The 
maximum likelihood principle implies that estimations of 
these parameters correspond to the maximum value of C 
(or, equivalently, ln£). The function Q can be rewritten in 
the form InC = -x /{2k) - (iVlnK)/2 - {lnw)/2 + const. 
When the weights Wi are fixed and known, the maximization 
of In yC can be performed elementary. The resulting value of 
6* is given by the least squares estimator. For the estimation 
of K, we obtain the well-known result k* = x^ {0*) / N . 

Let us now assume that we have r 'variance' parame- 
ters p (say, RV jitter of a given star observed with different 
instruments) entering in the model of Oi. Let us denote the 
full vector of d -I- r unknown (or at least poorly known) pa- 
rameters (0,p) as ^. For the square-additive model of cji, we 
should maximize In £ by and p simultaneously. The val- 
ues 6* and p* providing the maximum of In C, represent the 
joint maximum likelihood estimator (MLE) It is impor- 
tant that information about p* is automatically accounted 
for in the estimation 6* , and vice versa. An analytic max- 
imization of ln£ for the model ([l]) does not seem possible. 
However, an effective way of numerical maximization of the 
likelihood function will be described in Section [S] 

Any estimation is not of much use without associated 
uncertainty, i.e. without estimation of its variance. The 
variance-covariance matrix of an MLE is usually expressed 
using the Fisher's information matrix 

F(0=E(^^«^j=-E(^^j, (5) 

calculated for the true value of 4. The inverse represents 
an asympt otic (A'^ — » oo , i.e. large sample) approximation 
to Var C l|Lehman|[l98^ . § 6.4). In our case, the Fisher's 
information matrix can be written in the block form 

F=( f" l"")' (6) 

where the sizes of the submatrices match the dimensions of 
the vectors marked in subscripts. The calculation of F yields, 
in particular, that Fgp — and that fgg coincides with the 
Fisher's information matrix for the least squares estimator; 

Q^{fi'g(Sti'g/a^) (7) 

This implies that the vectors 9* and p* are asymptotically 
uncorrelated and the asymptotic variance-covariance matrix 
of 6* is the same as in the usual least squares approach. 



If our dataset is merged from several time series obtained 
at different observatories, we may be interested in separate 
estimations of RV jitter. It is not hard to show that these 
separate jitter estimations are asymptotically uncorrelated 
also. Finally, 

Var0* ~ Q'\ Varp* ~ e,^ = 2/{a-^)j, (8) 

where the index j means that the respective summation (*) 
should be restricted to the j"^ sub-dataset. The only seem- 
ing obstacle in practical use of (|5J comes from the fact that 
formally we should substitute the true values of parameters 

and p in these equations. In practice, we can substitute 
only the estimations 0* ,p* , which we have obtained before. 
This is admissible for calculating the asymptotic large sam- 
ple approximation of Var{*, because the estimations tend 
to the true values when TV grows. 

The MLEs possess many good statistical properties 
when the number of observations is large. Under certain 
regularity conditions, they are asymptotically (TV — » oo) 
unbiased (but see Section |6] for some cautio ns), asymptot- 
ically Gaussian and asymptotically efficient l|LehmarJll983l . 
chapter 6)0 The latter property means that the statistical 
uncertainties of MLEs approach the minimum possible ones 
when TV grows. Comparing the uncertainty of the MLE, pj, 
with the uncertainty of the MME from Section|31 we can ob- 
tain that their ratio is equal to ^ (cr^) (o-"*) /TV. Due to the 
Cauchy-Schwarz inequality, this quantity is not less than 1. 
This means that the MLE yields generally more accurate 
estimation of the RV jitter than the MME. For instance, 
the RV uncertainties of the Lick data for 51 Pegasi (see Sec- 
tion [TT]) imply roughly double advantage of the MLE. 

The MLE is organised so that the resulting value of 
X^/TV is always close to unity. This means that the x^ statis- 
tic can no longer be used as a measure of the fit quality. 
Instead, we should use some other statistic, based on the 
full likelihood function Q. For this statistic to be intu- 
itively clear, it should be measured in the same units as 
Vi. Therefore, it should be proportional to £~^^'^ . To find a 
suitable proportionality factor, let us assume for a moment 
that x^ = N exactly. Then ^-i/iv ^ o-gcome"'^ ^2^, where 
Cgoom is the geometric mean of tj,;. Therefore, already for 
the general case, we may introduce the following likelihood 
goodness-of-fit statistic: 

1 = r-'/'^e""- VV2^ ^ 0.2420£-'/'^. (9) 

This statistic describes naturally the overall scattering of RV 
measurements around a given RV model, for a given value 
of the RV jitter. 



^ This behaviour can be damaged in an incarefuUy chosen 
parametrization. For instance, it is a frequent case for hot Jupiter 
planets when the orbital eccentricity estimation looks like e = 
0.05 ±0.05 and the argument of the periastron uj is ill-determined. 
Then the distribution of e and u! is non-Gaussian. This is due to 
the formal singularity of the point e = in the polar coordinate 
system {e,u>). This trouble is easy to overcome by means of the 
change of variables x = ecosw, y = esinij. The joint distribution 
of {x, y) is already close to the bivariate Gaussian one. 
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6 BIAS REDUCTION 

It is well-known that linear least squares estimations are 
'unbiased', i.e. their mathematical expectations are equal 
to true values of parameters. This property is very impor- 
tant, because it allows us to hope that such estimations are 
related to true values at all. Both square-additive and multi- 
phcative models of RV uncertainties require non-linear like- 
lihood maximization to estimate the noise level parameter 
(k of aj). In general, maximum likelihood estima tions are 
biased, but their bias tends to zero as A*' — > oo (|Lehmanl 
1 19831 . § 6.4). Nevertheless, the biasing for real RV time se- 
ries with finite A'^ may become practically significant and 
may require a reduction. For instance, it is well-known that 
the maximum likelihood estimation k* = (0*) /N , derived 
in Section O is biased by 0{1/N) and the unbiased estima- 
tion is K* — {6*)/{N — d). In practice, we may quite have 
a set of d ~ 20 parameters of the Keplerian RV curve (for a 
four-planet system) with N ~ 100 observations only. In this 
case, the relative bias in k* (about d/N ~ 20%) exceeds the 
relative uncertainty of k* (about l/y/N ~ 10%). We may 
expect a similar biasing for cr J . The general reason of this bi- 
asing comes from the fact that residuals underestimate true 
errors in average. This underestimation increases when the 
number of free parameters grows. As an illustration, in the 
extremal case d — N we could plot a model curve transiting 
through all the data points exactly. In this case, all residuals 
would vanish. 

We need to reduce jitter bias so that the resulting esti- 
mation of would account for this reduction automatically. 
This reduction can be reached by means of proper C'(l/A'') 
modification of the functions (|4]) and This i s the ap- 
proach of 'preventive' bias reduction (iFirthlllQgsT ). For our 
specific goal, the likelihood function should be modified so 
that the residuals should be increased by the relative quan- 
tity ~ d/N, in order to reach a more accurate representa- 
tion of measurement errors. The following modification looks 
convenient in practice: 

ln£ = -xV(27) - (Ina) AflnV2^, (10) 

I = £-^/^e-0-VV2^ ~ 0.2420£-'/^, (11) 

where 7 = 1 — d/N. Clearly, such C'(l/A'^) modification 
should not destroy the large-sample properties (like asymp- 
totic normality and asymptotic efficiency) of the maximum 
likelihood estimator. But moderate- and small-sample prop- 
erties look now better. Maximizing (|10p instead of @ kills 
all bias in the estimation of k for the multiplicative model of 
RV uncertainties. In the case of the square-additive model, 
some residual bias may remain. This remaining bias can be 
calculated till the first order, 0{1/N), analytically, using 
cubic part of the Ta ylor expansion of ln£(^) near the true 
value of ^ (see, e.g., ICox fc Snelllll968l : lFirthlll993h . These 
calculations involve quite bulky tensor algebra and are omit- 
ted here. The final result (for the case r = 1) looks like 

(correction to crj) = i ^Tr (^Q"^q) - A-^^ , (12) 

where 

Q = <Me®Me/<^'), A = (a-2), = <0 . (13) 

The counterbalancing term in the equality (|12|l . containing 
d/N, was produced by our modification of the likelihood 



function. If all ai are equal to each other then the correc- 
tion (|12p is zero, as we could expect (recall that the mul- 
tiplicative model of ai is equivalent to the square-additive 
one for this case). When RV jitter is estimated separately for 
different components of the combined time series, we should 
apply this bias correction separately as well. In this case, 
it is necessary to restrict summations {*) in (|13p over the 
respective sub-datasets, but to keep the full summation for 
the matrix Q. The built-in correction provided by the like- 
lihood function modification (I10|l normally accounts for a 
large fraction of the bias in jitter. Therefore, the cross in- 
fluence of this bias on the estimations of is signiflcantly 
decreased. 

To c orrect the bias in estimations, the algorithm pro- 
posed bv lQuenouilld (|l956l ) may be used. This is also called 
the 'Jackknife' or 'leave-one-out' method. It is as follows; 

(i) Calculate the basic (biased by 0{1/N)) estimation x 
of a given parameter ^ from the full set of N observations. 

(ii) Construct N reduced time series with i*''(i = 
1,2, ...A'^) measurement omitted. Therefore, each reduced 
time series should consist of A' — 1 data points. 

(iii) Calculate A'' new estimations Xi{i — 1,2,... A'') of ^ 
by re-fitting with every of the reduced time series. The bias 
of x'i will be about ~ 1/(A'^— 1), hence these new estimations 
will be shifted with respect to x by about ~ (1/(A'^ ~ 1) ~ 
1/7V) = 0{l/N^). 

(iv) Calculate the sum bi = YliLii^'i~^)- The result bi = 
0{1/N) is the first-order bias of x. That is, the corrected 
estimation x — hi should be biased by C'(l/A'^) only. 

The main advantage of this algorithm is that its implemen- 
tation is model-independent and easy. Also, this algorithm 
does not require for the distribution of RV errors to be Gaus- 
sian. In addition, it can be directly applied to either 'vari- 
ance' (p) or usual (0) parameters. Unfortunately, it is rather 
time-consuming because it requires many non-linear fits. 



7 NON-GAUSSIAN ERRORS 

To write down the equality Q for the likelihood function, 
we have assumed that RV errors follow Gaussian distribu- 
tions. Some fears are sometimes expressed that RV errors 
in pl anet search surveys may be signi ficantly non-Gaussian 
fe.g.. iMarcv et al.ll2005l : iButler et al.|[2006 '). Then, strictly 
speaking, the function Q is not a likelihood function and 
estimations obtained from its maximization may be shifted 
with respect to the true MLE. Usually we have not enough 
information to construct the true likelihood function. Then 
the usage of simple Gaussian likelihood functions like Q or 
(|10|) may be reasonable. This is c alled somet imes the 'pseudo 
maximum likelihood' approach l|Bardlll974 § 4.18). 

How much the non-gaussianity of RV errors can affect 
the properties of the estimations obtained using the Gaus- 
sian likelihood function Q and its modification (|10p ? To get 
some preliminary answer to this question, let us consider a 
simplified situation of the least-squares algorithm from Sec- 
tion [3] with RV model being linear with respect to unknown 
parameters. The class of linear models incorporate, for in- 
stance, sinusoidal signals (C cos uit + S sin cut with a priori 
fixed frequency lj but free linear parameters C and S), and 
polynomial trends. This is the well-known linear regression 
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problem. The associated linear least-squares estimations can 
be expressed explicitly as certain linear combinations (or 
weighted sums) of the observations, regardless the shape of 
the input errors distribution. The general expressions for the 
coefficients are too unpleasant to be written down here, but 
they can be easily found in any textbook on the least-squares 
method. The errors of the derived linear estimations repre- 
sent just the same linear combination of the observational 
errors (again regardless the degree of their gaussianity) . This 
immediately implies the following properties of the linear 
least-squares estimations in the non-Gaussian situation: 

(i) If our RV model is correct, such estimations are ex- 
actly unbiased, regardless the shape of the distribution of 
the input RV errors. 

(ii) If the variances of the input RV errors exist (they may 
not exist, e.g., for heavy-tail Cauchy distribution) and are 
correctly modelled, the variances and correlations of derived 
estimations are exactly the same as in the case of Gaussian 
errors. In the non-Gaussian case, the linear least-squares 
estimator is no longer guaranteed to be strictly efficient, but 
still its variance is minimum possible among all unbiased 
linear estimators (the Gauss-Markov theorem). 

(iii) If the conditions of the central limit theorem for the 
given distribution of RV errors are satisfied, the joint distri- 
bution of the derived estimations tends to the multivariate 
Gaussian one when N ^ oo. 

The mentioned general properties of the le ast-squares es- 
timat ors are well-known in statistics (e.g., iKoroluk et al.l 
1 19781 . §23.2.6). 

Of course, the models of the RV curve met in planet 
search syrveys typically incorporate non-linear Keplerian 
RV functions. We should not expect that the nice prop- 
erties of the linear least-squares estimations with non- 
Gaussian input errors should hold true for the more com- 
plicated non-linear pseudo maximum likelihood case. How- 
ever, we can suspect that at least some of these properties 
may be conserved approximately in the asymptotic sense 
for — > oo. This problem was consi dered rigorously by 
iGourieroux. Monfort fc TrognonI (|l984l ) (pay particular at- 
tention to their Section 6). One may be surprised, that (of 
course under certain regularity conditions) many important 
asymptotic properties of maximum likelihood estimators, 
constructed for Gaussian errors, are conserved in the pseudo 
maximum likelihood case, i.e. when the errors do not follow 
Gaussian distributions. For example, the pseudo maximum 
likelihood estimators are asymptotically unbiased and Gaus- 
sian. However, the asymptotic efficiency may be lost: we 
cannot construct even asymptotically efficient estimator if 
the shape of the distributions of the RV errors is not known 
precisely. The asymptotic variance-covariance matrix of es- 
timations 0* ,p* in the case of non-Gaussian errors can be 
derive d from the formulae given in the Append ix 5 of the pa- 
per by IGourieroux. Monfort fc TrognonI l| 19841 '). The matrix 
Varfl* ~ is unchanged (in the asymptotic large-sample 
approximation). Jitter estimations corresponding to differ- 
ent observatories are uncorrelated again. However, a non- 
zero skewness of RV errors inspires some correlation between 
6* and p*, and an excess kurtosis distorts the variances of 



pV- 

Varp* ^ ,^ + 1^1^^^. (14) 

It is important that if a large skewness (i.e., asymmetry) 
of RV errors would be checked to be negligible, large cross- 
correlation between 6* and p* should not be expected. The 
variances of Pj may either increase (for leptokurtic RV er- 
rors. Ex > 0) or decrease (for platykurtic RV errors. Ex < 0). 
Note that the expressions (|14|) do not require for the shape 
of distribution of the RV errors to be known in advance. 
Provided only the skewness and kurtosis are known, it is 
possible to use these expressions in practice. For instance, 
if the kurtosis of RV errors is constant, the variance of the 
corresponding jitter estimation should increase by the fac- 
tor (1 + Ex/2). The expressions (|14|) were checked by means 
of Monte Carlo simulations, assuming different simple non- 
Gaussian distributions for simulated RV errors (e.g. uniform 
one). The predictions of analytic formulae (|14p for the jitter 
estimation variance were found to be in an excellent agree- 
ment with results of numerical simulations, at least for 
as big as a few hundred. 

Of course, we should satisfy certain conditions of regu- 
larity for the theoretical results described above to be appli- 
able. The rigorous formulation of these conditions is given i n 
the Appendix 1 bv IGourieroux. Monfort fc TrognonI l|l984h . 
These incorporate: 

(i) Certain requirements of boundedness and integrability 
for the distribution of the RV errors (roughly speaking, too 
heavy tails are not allowed). 

(ii) Conditions of smoothness and boundedness for the 
equations of the model (RV) curve (and also for the model of 
variances, but our square-additive and multiplicative models 
of uncertainties are very simple and certainly satisfy them). 

(iii) Requirement that (roughly) any single observation 
should not strongly affect the final estimations. This put 
certain condition of 'naturality' on the sequence of observa- 
tional timings and statistical weights (and also on the RV 
model) . 

In fact, these conditions are not qualitatively new. The first 
one originates from the requirement of the central limit teo- 
rem. The second one originates from the regularity condi- 
tions from the maximum likelihood estimations theory. The 
third one originates from both fields. 

However, I could not find enough information in the 
literature about skewness and kurtosis of RV errors in cur- 
rent planet search surveys. By this reason, I could apply only 
formulae (jSj, valid for Gaussian RV errors, when calculating 
the uncertainties of estimations in Section 1111 It is impor- 
tant to note that still the degree of possible non- Gaussianity 
of the RV errors in planet searches is not cl early estimated. 
In a r ecent study of the Keck RV survey, ICumming et al.l 
(■2008') did not reveal clearly any strong non-Gaussianity. 

Non- Gaussian errors may lead to extra C'(l/A'') biasing 
of 0* and p* . This extra bias can be calculated till the first 
order using the same approach based on the Taylor expan- 
sion of ln£({) as in Section [6] A non-zero skewness of RV 
errors leads to an extra bias ~ As/N in estimations of G 
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(including estimations of planetary parameters). Consider- 
ing that very large skewness of RV errors is unlikely, the 
'Gaussian' part of the bias in should dominate in practice. 
Therefore, the bias in G inspired by non-Gaussain errors 
should not be a practical trouble. A kurtosis excess of RV 
errors adds some extra bias ~ Ex/N in the estimations of 
jitter. However, this effect is expected to be negligible even 
for the kurtosis excess as large as Ex =1 — 3. For instance, 
this bias vanishes when the kurtosis is constant. In any case, 
the first-order bias in parameters G can be removed by the 
Quenouille's algorithm (see Section [6]). 



8 NUMERICAL CALCULATION 

It is necessary to propose numerical algorithms performing 
maximization of the function (|1U|I . For the sake of simplicity, 
let us put r — 1, p = a^. The extension to the case r > 1 
will be straightforward and easy. Let us consider the function 
g = const —2 In £ to be minimized by G and p: 



In 1 + ■ 



+ 



7('7mcas,i + P) 



(15) 



' mcas , i • 



This function is formally defined for p > po = — min cr„ 
Note that negative values of p are not senseless. They indi- 
cate that the instrumental uncertainties specified are in fact 
overestimated. 

It seems better to minimize g{G,p) in two steps. In 
the first step, we will obtain a second-level target function 
h{G) = mirip g{G,p) — g{G,p*{0)), where p*(0) denotes the 
value of p for which this minimum is achieved. It is evident 
that in non-degenerated situations \imp ^ + oo g{G ,p) — +oo 
and liriip pg g{0,p) = +oo. Therefore, for any fixed G at 
least one minimum by p exists. The one-dimensional min- 
imization by p for a fixed can be precisely and rapidly 
performed by simple Newtonian-like algorithms. 

A robust situation with only one solution p* (0) > po for 
a given usually takes place. However, sometimes we may 
deal with the following ill conditioned case. Suppose that for 
some — 00 the residual corresponding to (Jmoas.i ~ —po 
vanishes. Then the function (|15p has no global minimum. In 
this case, \imp-,pg g{Go,p) = — oo and the respective solu- 
tion G — Go and p = po is not physically sensible. For well- 
conditioned cases, real minimization algorithms converge to 
good solutions which are far from these singularities. The 
situations when a numerical algorithm falls in the singular- 
ity are very seldom in practice and appear when the RV 
uncertainties span a wide range and/or the RV models are 
overloaded (contain too many free parameters) and/or they 
are close to being degenerate. These cases represent a nu- 
merical problem and should be identified during the mini- 
mization. A simple test 1 -I- p/cr^oas,i < 0.01 is sufficient to 
diagnose almost all the singular cases. 

In the second step, the function h{G) should be mini- 
mized. This may be performed by standard non-linear least 
squa res algorithms like the Levenberg-Marquardt-Gauss one 
ifSard 1974,, §§ 5.8-5.11). To show this, we need to check that 
the gradient and the Hessian matrix both can be calculated 
in the same way as during the minimization. Firstly, the 
gradient h'{G) is equal to the partial derivative g'g{0,p* (0)). 
Within the factor 7, the last partial derivative is the usual 



gradient of the function (calculated for the jitter p*(0)). 
Secondly, we need to check that the Hessian matrix h"{G) 
can be calculated using the Gauss' approach. Recall that the 
full Hessian matrix for the function X^(^) is given by 



(16) 



The second term in (fT6|l has magnitude 0{VN) and is ne- 
glected in comparison with the first one, which has the mag- 
nitude 0{N). This is the commonly used Gauss' approach 
which allows the calculations of second-order derivatives of 
H to be avoided. It can be shown that the same approxima- 
tion is valid for the matrix 7/1" (0). The exact expression of 
'yh"{0) contains extra terms having magnitude 0(1) (i.e., 
0(A''")) only, that is even less than the second term in (|16|) . 



9 TESTING HYPOTHESES 

Often we need to choose between at least two hypotheses, a 
base one TL and an alternative one tC, based on the RV data. 
Usually these hypotheses are defined by some parametric 
temporal models of the RV curve, /iw(t, 0'h) and yLKit, Gk)- 
Here vectors Gn and 0)c contain dn and djc unknown pa- 
rameters. We will assume that Ti is nested in K., that is 
0K = {Gh,0} and iiK{t,GK.) = ^l■H{t,G■H) + fi{t,0) where 
d = d)c dn quantities parametrize the model /i of some 
extra RV variability. The parameters G are chosen so that 
this extra signal vanishes when = 0: /i(t, = 0) = 0. We 
wish to test whether the hypothesis TC : G — (no signal) 
is consistent with our RV data or it should be rejected in 
favour of the alternative K, : G ^ Q (signal exists). The pa- 
rameters 0k: are supposed to belong to some domain Qk. in 
d]c dimensions. The condition = cuts in this domain a 
hypersurface Qn of dimension dn < die- Thus we can refor- 
mulate our goal as to check, whether the hypothesis £ Ot^ 
is consistent with the RV data or it should be rejected in 
favour of the alternative £ Ok; \ Qn ■ 

There are many practical tasks which can be embedded 
in this mathematical framework. For instance, often we need 
to test existence of an extra periodic RV variation of a given 
frequency or an extra long-term RV trend. The possible ex- 
tra periodicity may be modelled as a sinusoidal harmonic, 
and the possible trend as a linear or quadratic function. 

The common tools used to solve such problems are the 
and F tests. The x^ test is based of the difference between 
the x^ functions, calculated for the best-fitting RV models 
for the hypotheses H and IC. To apply the x^ test, we need 
to know the full RV uncertainties cr;. The F test is based 
on the ratio of the same x^ functions. The F test is more 
fiexible than the x^ one: it can process cases when only the 
weights Wi are known a priori, and the RV uncertainties are 
calculated according to the multiplicative model ([2]). The 
factor K is estimated implicitly in the F test. In our case, the 
RV uncertainties are given by the square-additive model ((T|, 
and the F test cannot be applied. The RV jitter cr* has 
to be estimated explicitly. Doing so, we can construct the 
logarithm of the likelihood ratio statistic 



max In £. 



max ln£j| 



(17) 



Here, the maximization of In £ by p means that the RV jitter 
are estimated explicitly, together with the usual parameters 
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of the RV curve, 0. The resulting best-fitting values are then 
used to construct the logarithm of the ratio of the maximized 
likelihood functions corresponding to hypotheses TC and IC. 
For the purposes of bias reduction, it is better to use the 
following modification of the likelihood ratio; 



max In Ck — max In Cn 



where Nn = N — dn and A'^^ — N — d/c- The modified 
likelihood functions ln£ are different for hypotheses Ti and 
IC, because they contain different correctors y-n = Nn/N 
and "fic = N/c/N. Note that if the multiplicative model 
were assumed for ai, the functi on (1181) would coincide with 
the statistic Z3 from the paper (|Baluevl|2008al ) . The square- 
additive model of at generates another form of Z, which 
is preferred for testing statistical hypotheses in RV planet 
search surveys. Note that definitions (|17|l and (|18|> do not 
require strict linearity of the models. 

A large value of the statistic Z indicates that the base 
hypothesis may be wrong and the specified alternative model 
is more realistic. However, random RV errors may also pro- 
duce similar values of Z. To compute statistical significance 
of the observed value of Z we should know the distribu- 
tion of Z under the base hypothesis Ti. Of course, there is a 
little hope that this distribution can be calculated exactly. 
Nevertheless, many asymptotic (A*' 00) results are known 
for the likelihood ratio statistic. In particular, the distribu- 
tion of the quantity 2Z (as well as 2Z) converges to the 
distribution with d degrees of f reedom, if certain regular- 
ity c onditions are satisfied (e.g. iProtassov et al.ll2002l : ISenI 
Some of these regularity conditions are technical and 
are satisfied in the majority of applications. However, other 
conditions may not be satisfied in many practical cases, and 
therefore they deserve to be checked before applying the 
asymptotic distribution to Z. It is worth noting that: 

(i) The spaces of parameters should be nested. Oh C Qic- 
This requirement is already built in our formulation of the 
hypothesis testing problem. 

(ii) The subspace Q-h should lie in the interior of Ojc. It 
should not lie on the boundary of Ok. Otherwise, the asymp- 
totic distribution of the likelihood ratio st atistic is not the 

d istribution with d d egrees of freedom |Protassov et al.l 
|2002| ). See the paper bv ISelf fc Li"an3 (Il987l ) for a general 
algorithm of constructing the asymptotic distribution of Z 
(or Z) in this non-standard case. Typically, when Oh he 
on the boundary of O/c, the asymptotic distribution of the 
likelihood ratio statistic appears to be some mixture of x^ 
distributions with different numbers of degrees of freedom, 
but more complicated cases are also possible. 

(iii) Equations of the RV models, ij.n(t,0n) and ^{t,6), 
should satisfy certain conditions of smoothness and bound- 
edness. 

Note that the same (or similar) regularity conditions are 
equally required to hold true when using the F test and the 
X^ test as well. These conditions may not to hold true for a 
given parametrization but simultaneously may be satisfied 
for some other one. The likelihood ratio statistic and its dis- 
tribution are invariant with respect to a re-parametrization. 
Hence, it is sufficient for the regularity conditions to hold 
true for only one parametrization. 

Suppose that all the necessary conditions are satisfied. 



and the distribution of 2Z indeed converges to the x^ one. 
This convergence is not uniform. Larger values of Z corre- 
spond to larger displacements in the parameter space. These 
increase non-linear effects and require larger A''. I have fol- 
lowed the convergence of Z to the x^ distribution for the 
square-additive model ((T| by means of Monte-Carlo simu- 
lations for various structures of time series and simple RV 
models. The simulations yielded the following empirical con- 
vergence condition: 



N > pZs (valid for A^ > 30, FAP > 10"^ s < 1), 



(19) 



where s — ln(max cr/ min a), p = 2.5 — 3 when testing RV 
model 'constant velocity' vs. 'constant -I- linear trend' (dn = 
d = 1), and p = 5 — 6 for RV models 'constant' vs. 'constant 
+ sinusoidal harmonic with a fixed frequency' {dn = l,d = 
2). The condition (|19|) is not stringent. The extra RV jitter 
softens differences between ai, so that s < 1 usually. Note 
that we are practically interested in the values Z = 3 — 10 
producing FAP = 10"^ - 0.1 (for small d). 



10 LIKELIHOOD PERIODOGRAMS 



The iLombl l|l976l ') 4Scargi3 l| 19821 ') periodogram and its nor- 
malizations require fixed weights of observations; hence, 
they assume multiplicative model for RV uncertainties. For 
planet searches, it is preferred to use a periodogram with 
built-in estimation of the square-additive RV jitter. Such pe- 
riodogram may be constructed from the modified likelihood 
ratio s tatistic Z in the same way as it is described in (|Baluevl 
l2008al ) for the statistic. The resulting periodogram Z{f) 
will be a function of the frequency / of a trial periodic RV 
signal (modelled by a sinusoidal harmonic or by some other 
periodic function), so that every single value Z{ f) represents 
the statistic (|18|) . 

To assess statistical significance of candidate periodici- 
ties, we should know the distributions of the likelihood ratio 
periodograms. In practice, orbital period of a possible planet 
is not known a priori and a wide frequency range is scanned 
in order to find the maximum periodogram peak. By analogy 
with single- value distributions, we can expect that distribu- 
tions of maximum values of Z{f) converges (for A'' — > 00) to 
the distribution of the maximum of the least squares peri- 
odogram z{f) from (|Baluevll2008"al '). Here we should be care- 
ful, because the convergence of distributions of maximum 
values of periodograms is not necessary achieved in practice 
(|Schwarzenberg-Czernvl [19981 ) . Recall that the conv ergence 
condit ion for the periodogram 2:3 (/) from the work (iBaluevI 
l2008al ) is 23 <^ N. This periodogram is a special case of the 
likelihood ratio periodogram for the multiplicative model of 
RV uncertainties. Therefore, we can expect that our like- 
lihood ratio periodogram with a built-in jitter estimation 
requires a convergence condition similar to (|19p . Results of 
Monte-Carlo simulations are consistent with this assump- 
tion if p = 5 — 6 or even less (depending on the aliasing 
and on the desirable precision of false alarm probability). 
Unfortunately, this problem appears too complicated to be 
discussed here in more detail and probably requires a sepa- 
rate investigation in future. 
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11 APPLICATION TO REAL DATA 



Table 1. RV jitter for several stars with planetary systems. 



The mathematical tools described above may be applied 
to real RV data from planet search surveys. For this pur- 
pose, I have selected several planetary systems with well- 
determined orbital configurations and large series of observa- 
tions from d ifferent observatorie s. RV data were taken from 
the works bv lButler et~all (120061 ) for t he stars 51 Peg, 70 Vir, 
14 Her, HD83443, 54 Psc, /x Ara; by iNaef et al.l (|2004|) for 
51 Peg, 70 Vir, 14 Her, 55 Cnc: by iMavor et all (|2003 ) 



for HD8344 3: bv IWittenmver. Endl fc Cochran 



14 He r; by iPepe et al 



1I2OO6I) for HD6983q ; by 



(120071) for M Ara; by 



il2007F 



for 



Lovis et aU 



Wittenmveret al.' ('2007|) f or 54 Psc 



bv iMcArthur et all (|2004l ) and bv iFischer et all l|2008h for 
55 Cnc. 

For almost every of these stars we have two or more in- 
dependent RV datasets. The configurations of orbits in these 
systems are determined reliably. For each star, we can write 
down the model of the RV curve containing a number of free 
parameters to be estimated. The joint estimation of the RV 
curve parameters and of the RV jitter was performed using 
maximization of the modified likelihood function (|10|l . The 
resulting jitter estimations are shown in the fourth column 
of Table [1] These estimations include analytic bias correc- 
tion (|12p . though this correction often appeared negligible 
due to a large number of observations. The corresponding 
estimations of planetary parameters are omitted (perhaps, 
they are not so interesting here). We can clearly see that 
RV jitter for distinct instruments may differ largely, even 
for one and the same star. Typically, RV jitter for ELODIE 
and CORALIE spectrograph's are significantly higher than 
for other instruments. The only exception is the star 70 Vir, 
for which the ELODIE RV jitter is consistent with zero (two- 
sigma upper limit is 3.2 m/s). In contrast, the HARPS jitter 
are remarkably smaller, dropping sometimes below 1 m/s 
level. Data from Keck, Lick and Angl o-Australian o b serva - 
tories demonstrate intermediate cases. iFischer et al.l (|2008l ) 
estimated RV jitter for Lick and Keck data for 55 Cnc by 

3.0 m/s and 1.5 m/s. However, the corresponding values 
from Table [T] are 5.19 m/s and 4.33 m/s. This indicates 
some extra RV instability having standard deviation about 

4.1 m/s. The actual nature of this extra RV instability is un- 
clear. Surprisingly, the RV jitter for the HJS data for 14 Her 



is definitely negative (ct^ 



-(4.5 m/s) ). This indicates 



that the RV uncertainties of these me asurements, quoted by 
IWittenmver. Endl fc CochranI (|2007l ). are overestimated by 
about 20%. More careful investigation reveals that the r.m.s. 
of the RV residuals for the HJS data of 14 Her is 5.85 m/s. 
This value of r.m.s. is quite reliable, because the RV model 
is well constrained by the data from two independent teams 
(from Lick and ELODIE) . Simultaneously, the stated instru- 
mental RV uncertainties of HJS are ranged between 6.4 m/s 
and 12.7 m/s with the geometric mean of about 7.4 m/s. 

Large RV jitter for ELODIE and CORALIE data in- 
dicate that some systematic errors are not accounted for 
in their RV uncertainties. Partially, these systematic errors 
can be explained by extra annual RV variations. This can 
be established by means of including an extra harmonic 
term ^cos(27r(t-r)/lyr) = Ccos(27rt/lyr)-|-S'sin(27rt/lyr) 
in the RV model. Here, the parameters C, S (equivalently 
A, t) were estimated from the RV data together with plan- 
etary parameters. Then it can be checked, whether the 



y Lciii 


Instr.* A'^ nomin. a^,^ 






51 Peg3 
b 


ELD 153 9.93(88) 
LCK 256 0.53(2.37) 


11.2(1.2)Jun 28.5(9.3) 


6.56(83) 
0.45(2.81) 


70 Vir 

b 


ELD 35 
LCK 74 


-3.0(1.6)2 
4.22(80) 


5.1(1.6) May 16(16) 


-2.8(1.8)2 
3.58(79) 


14 Her^ 
b 


ELD 119 
KCK 49 
HJS 35 


7.32(96) 
1.25(44) 
-4.46(84)2 


6.4(1.6) Jan 6(12) 


6.45(96) 
1.23(45) 
-4.44(86)2 


HD69830 
b, c, d 


HRP 74 


0.22(26) 






HD834433 
b 


COR 257 
AAT 23 
KCK 28 


6.83(49) 
5.9(1.7) 
2.74(57) 


5.8(1.4) Oct 12.5(9.5) 


6.40(48) 
6.0(1.7) 
2.70(58) 


54 Psc 
b 


LCK 121 
KCK 42 
HET 83 


5.78(52) 
3.86(52) 
10.1(1.3) 


4.73(90) May 27(13) 
11.7(3.1) Dec 18(27) 


5.87(54) 
2.30(40) 
7.8(1.0) 


/I Ara 
b, c, d, e 


COR 40 
AAT 108 
HRP 86 


5.67(95) 
2.37(25) 
1.52(15) 






55 Cnc5 

b, c, d, 
e, f 


ELD 48 
LCK 250 
KCK 70 
HET 119 


14.2(1.9) 
5.19(33) 
4.33(36) 
6.95(86) 


8.3(3.5) Oct 21(34) 
8.2(1.8) Mar 14(15) 


13.3(1.8) 
5.19(33) 
4.33(37) 
5.22(74) 



The uncertainties of the estimations are given in the parentheses 
in the units of last digits. For instance, 7.32(96) means 7.32±0.96, 
and 3.0(1.6) means 3.0±1.6. The quantities A and r are the semi- 
amplitudes and the epochs of the maximum RV of the best-fitting 
sinusoidal annual drifts. 

^ Uncertainties of jitter estimations are calculated in the linear 
Gaussian approximation as 5 = £/(2(T*), where e is the asymp- 
totic uncertainty of the estimation of trj . When CTj- is comparable 
with (or less than) (5, its distribution is far from Gaussian. Then 
it is necessary to return back to the quantity cr^ (not affected by 
the degeneracy) and to its uncertainty e = 2(5cr* . 
2 Negative values of symbolically reflect that corresponding 
estimations of crj are negative. 

^ A linear/quadratic trend was included in the RV model. 

ELD = ELODIE, COR = CORALIE, HRP = HARPS, LCK 
= Lick Observatory, KCK = Keck Telescope, AAT = Anglo- 
Australian Telescope, HET = Hobby-Eberly Telescope, HJS = 
Harlan J. Smith Telescope. 

^ Dynamical fit, taking planetary perturbations into account. The 
orbits were assumed coplanar and seen edge-on (other moderate 
inclinations gave slightly difi'erent but similar results). 



maximum value likelihood function is increased statisti- 
cally significantly, that is whether the observed value of the 
corresponding statistic Z is statistically significant. As we 
have d — 2 (two free parameters of the annual harmonic), 
the asymptotic distribution of Z (and Z) is exponential: 
FAP = Pr{^ > z} « e"^ Columns 5,6,7 in Table [l] con- 
tain the best-fitting parameters of the annual variations and 
residual RV jitter, for the cases when the significance of an- 
nual terms was large enough. For the star 51 Pegasi, for 
instance, the value of the statistic Z (almost equal to Z in 
this case, due to a very large N = 409), associated with 
the annual periodicity, was about 35. For comparison, the 
3-sigma significance level is Z ~ 5.9. 
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Figure 1. Top: likelihood residual periodogram of ELODIE ra- 
dial velocities of 51 Pegasi. The RV oscillation induced by the 
planet 51 Peg b was included in the base model fi-}^. For the 
model fj, of the extra signal to be tested, the harmonic function 
was adopted (as it is done in the Lomb-Scargle periodogram). The 
clear peak with Z = 30.8 at the period of 359 days most likely 
corresponds to annual instrumental errors. Three strong peaks 
near one day (Z = 27.3 at 23''52"'6'', Z = 23.7 at 24^0^5^, and 
Z = 18.4 at 24'^4'"2^) may be interpreted as aliases. Bottom: 
the periodogram with base RV model 'planet b + linear trend + 
annual periodicity'. This periodogram is clean. 



It is interesting to consider planet candidates having 
orbital periods consistent with one year. There are about 
fifteen such planets. The masses of all of them but one are 
rather large (several Jupiter masses), producing large RV 
semi-amplitudes ~ 100 m/s. However, the planet HD74156 
d announced rec ently on the basis of HET observations 
d Sean et al.ll2008l ) has relatively small mass and induces the 
RV oscillation of ~ 10 m/s only. Given similar amplitude of 
the annual periodicity in the HET data for the stars 55 Cnc 
and 54 Psc, the discovery of HD74156 d looks suspicious. 
More careful analysis shows that RV data for HD74156 from 
ELODIE ( Naef ct al. 2004) also show an annual variation of 
20 m/s, but in opposite phase. This inspires strong doubts 
about the existence of the planet HD74156 d. This planet 
candidate could quite be a false detection made due to an- 
nual errors in RV data from HET. In any case, its orbital 
parameters may be strongly dis torted and are unreliable. 

The work (|Baluevl [2OT8bl ) provides a more intricate 
example, dealing with RV data for the system around 
HD37124. It presents a careful application of the method- 
ology described here. It further demonstrates either the im- 
portance of the differences between the effective RV jitter 
for ELODIE, CORALIE, and Keck data, or the importance 
of annual errors in obtaining suitable orbital solutions. 



We can see that annual RV errors in planet search 
surveys represent a frequent phenomenon (although cases 
free from these systematic errors are also not rare). Pos- 
sible sources of such systematic errors may have various 
nature. For instance, for ELODIE and CORALIE they 
may be partially inspired by weak telluric lines which were 
not completely excluded in old cross-correlation templates 
(M. Mayor, private communication). Keck and Lick data 
published before 2006 contain an inaccuracy due to non- 
relativistic barycentric correction. This inaccuracy was be- 
low 1 m/s for the majority of stars, but in rare cases reached 
3 — 5 m/s (J.T. Wright, private communication). Sometimes 
annual periodicities can produce rather strong peaks on pe- 
riodograms (Fig. [T}. However, a more dangerous situation 
takes place when such contaminating RV variations are not 
seen on the periodogram clearly. Despite of this fact, they 
can produce significant distortions of the best-fitting orbital 
model of the planetary system. Therefore, often it may be 
useful to check how much the best-fitting orbital configura- 
tion of a planetary system is changed if the annual harmonic 
ter m is added to the RV model. 

iButler et al.l (|2006| ) report on the detection of an extra 
RV trend of -1.64 ± 0.16 m/(s-yr) in the Lick RV data for 
51 Peg. This trend may be induced by an extra unseen com- 
panion on a long-period orbit. However, the ELODIE data 
alone yield an estimation of —0.15 ± 0.40 m/ (s-yr), which is 
poorly consistent with the Lick data. Only after addition of 
an annual harmonic to the model of the ELODIE data the 
long-term trend can be confirmed. Significance of this slope, 
based on the ELODIE data only, now corresponds to 2.7- 
sigma level. Its magnitude of —0.94 ± 0.34 m/ (s-yr) is now 
much better consistent with the estimation bv IButler et al.l 
l|2006t ). although some residual difference at the level of less 
than two sigma may indicate some extra (probably non- 
periodic) systematic RV errors, yet to be corrected or taken 
into account. The joint estimation, based on the Lick and 
ELODIE data, yields -1.46 ± 0.16 m/(s-yr). 



12 CONCLUSIONS 

In this paper, the problem of poorly constrained radial ve- 
locity jitter in planet search surveys is considered. An al- 
gorithm for the RV curve fitting with a built-in accounting 
for this jitter is developed. In many cases, this algorithm 
gives much better accuracies of estimation of the full RV 
jitter than the usage of only astrophysical jitter estimations 
based on the empirical correlations with stellar character- 
istics. This algorithm is based on the maximum likelihood 
principle and includes a series of bias corrections, either an- 
alytic or numerical. An extension of the Lomb-Scargle pe- 
riodogram with a built-in jitter estimation is proposed. An 
effect of non-Gaussian RV errors is discussed and is shown to 
be tolerable for large time series. It was shown that numeri- 
cal computations based on this method may be implemented 
by means of the usual least squares Levenberg-Marquardt- 
Gauss algorithm with slight modifications. 

This methodology can be useful for obtaiinng best- 
fitting orbital configurations in extrasolar planetary sys- 
tems, especially in the case when inhomogeneous RV data 
(taken from different observatories) are merged in the anal- 
ysis. It can also provide better input data for improving the 
empirical correlations found in activity models, since the 
maximum likelihood approach provide a better accuracy of 
jitter estimations than the tradiationally used method of 
moments (which is based on the difference between observed 
scatter of RV residuals and instrumental noise level). 

The application of these mathematical tools to several 
stars with published RV data revealed that the effective val- 
ues of the RV jitter may be quite different for different in- 
struments (even for one and the same star). Probably, the 
main reason of these differences is a poor knowledge of in- 
strumental RV errors. In many cases (especially for ELODIE 
and CORALIE data), an extra annual harmonic in the RV 
model leads to significant improvements in fit quality and 
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somewhat decreases too large cfloctive RV jitter. It is shown 
that the planet candidate HD74156 d with orbital period 
close to one year may be a false detection made due to an- 
nual errors in the HET RV data. 
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APPENDIX A: NOTATION 

Let us introduce the following operation: 

JV 

(0(t)) = ^<^(ii), 

i=l 

where ti is an i"^ observational epoch. Similar summation 
(ffli) may be defined for a discrete sequence at. For shortness, 
the argument t and the index i are omitted in the text. 

All vectors (including gradients) are assumed to be col- 
umn ones by default. The notation {xi, X2, ■ ■ ■} correspond 
to a vector constituted by elements of the vectors xi,X2, ■ ■ ■■ 

If X, y are vectors then x ®y := xy^ is a matrix con- 
stituted by the pairwise products Xiijj. 

Var X :— K{x ® a;) — Ea; ® Ea; is the variancc-covariance 
matrix of the random vector x, and Gov(x, y) := E(x(8)y) — 
Ea; (8) Ey is the cross-covariance matrix of x and y. 

This paper has been typeset from a T^jX/ file prepared 

by the author. 



